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SUMMARY 

A  numerical  scheme  is  outlined  for  the  calculation  of  two  dimensional  time-linearised 
transonic  flow  about  an  aerofoil  or  collection  of  components,  using  the  Green's  function 
method.  Computed  results  are  presented  for  various  configurations  for  the  simpler  sub¬ 
problem  of  steady  subsonic  flow. 
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1.  INTRODUCTION 


This  article  presents  tin  outline  of  a  numerical  scheme  for  calculating  two- 
dimensional  unsteady  transonic  flow  about  an  aerofoil,  or  a  collection  of  aerofoil 
components.  The  aerofoil  components  are  assumed  to  be  thin,  and  the  effects 
of  viscosity  are  neglected.  Hence,  the  flow  is  assumed  to  be  governed  by  the 
transonic  smtdl  disturbance  equation  [1],  [7], 

The  approach  followed  here  is  the  Green’s  function  (integral  equation)  method. 
This  method,  see  [8-10]  and  also  [2],  was  first  introduced  to  compute  unsteady 
subsonic  aerodynamics  in  three  dimensions.  However,  the  method  was  ex¬ 
tended  in  [14]  to  also  handle  transonic  aerodynamics  (see  also  [3]).  The  main 
advantage  of  the  method  is  that  it  is  easy  to  implement  for  complex  geometries, 
since  body-fitting  grid  systems  are  not  required  [5], [15]. 

In  [14]  and  [3],  the  Green's  function  method  involves  “time- integration”.  That 
is,  the  flow  about  the  aerofoil  is  ctJculated  at  each  moment  as  the  program 
steps  forward  in  time.  Alternatively,  significant  savings  in  computer  time  can 
be  made  by  “time-linearising”,  i.e.  by  assuming  that  the  aerofoil  oscillates 
harmonically  about  some  metin  steady  position  ([5],  [8],  [9],  and  [11]). 

The  cost  of  “time-Unearisation”  is  an  inability  to  model  highly  non-linear  phe¬ 
nomena  such  as  shock  appearance  and  disappearance,  and  large  scale  shock 
motion,  such  as  that  identified  in  [13].  However,  the  method  is  compatible 
with  conventional  subsonic  (linear)  flutter  analysis,  and  significantly  easier  to 
program  than  “time-integration”  methods. 

As  a  first  step  towards  a  full  “time-linearised”  transonic  capability,  results  are 
presented  here  for  a  simple  sub-problem  of  two-dimensional  unsteady  tremsonic 
flow,  namely  the  problem  of  two-dimensional  steady  subsonic  flow.  This  pre¬ 
liminary  problem,  for  which  analytical  solutions  tire  available,  provides  a  useful 
initial  check  on  the  method’s  accuracy  and  a  pointer  to  possible  problems  in 
the  full  transonic  formulation. 


2.  TRANSONIC  PROBLEM 

Under  the  assumptions  that  the  fluid  is  isentropic  and  inviscid,  the  flow  is 
irrotational  and  the  aerofoil  is  thin  and  undergoes  small  unsteady  motion, 
there  exists  a  perturbation  potential  (t>  which  must  satisfy  the  unsteady  small 
perturbation  trtinsonic  equation  [1],  [7],  namely 

[(1  -  M^)  -  A/2(1  +  7)^xj«^ix  +  <l>yy  -  -  M^4>U  =  0.  (1) 

where  7  is  the  ratio  of  specific  heats  and  the  x-axis  is  in  the  streamwise  direc¬ 
tion.  The  Mach  number  M  is  assumed  to  be  nem,  but  less  than  unity,  meaning 
that  only  “sub”-transonic  flow  will  be  considered  here. 

In  order  to  reduce  the  computational  time  and  complexity  of  the  formulation, 
the  governing  equation  is  “time-linearised”  as  in  [5]  and  [11].  That  is,  the 
potential  is  separated  into  a  steady  part  and  a  small  harmonically  oscillating 
unsteady  part.  This  leads  to  a  non-linear  boundary  value  problem  for  the 
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steady  potential  and  a  lincM  boundary  value  problem  for  unsteady  potential 
(which  is  linked  to  the  steady  problem). 

Equation  (1)  is  “time- linearised”  by  setting  4>  =  4>s  +  ‘t>u,  where  «iu  C 
and  keeping  only  the  higher  order  terms.  Accordingly,  the  steady  potential  4>s 
must  satisfy 

(1  —  M‘^)(I>szx  +  <i>syy  —  +  'y)ipsx4>sxx  (2) 

which  can  be  rewritten,  using  the  Prandtl-Glauert  trtinsformation  A'  =  x, 
Y  =  /3y  (setting  3  =  \/l  —  M'^),  as 


<i>sXX  +  ‘I’sYY  ^ 

where  the  non-linettr  term  <7*  is  given  by 

M2(H-7),  , 

“  1  -M^  <PsX<t>sXX- 


(3) 

(4) 


The  unsteady  potential  (pu  must  then  satisfy  the  linear  equation 

(1  —  M^)0uxx  +  4‘uyy  ~  2A/^(^uij  —  M  4>utt  ~  -f  7)^^[<^ui'^sxl.  (5) 

Application  of  the  above  Prtindtl-Glauert  transformation,  plus  the  further  sub¬ 
stitution  (assuming  harmonic  time  dependence)  of 

=  (6) 


into  (5)  leads  to  an  equation  for  the  new  dependent  variable  $u,  namely 

^uXX  +  ^uvy  +  =  o-u,  (7) 

where  h  =  kM/0  and  the  <7u  term  on  the  right  hand  side  (which  is  linear  but 
dependent  on  the  steady  potential)  is  given  by 


<7u 


M2(1-|-7)  d 
1-Af2  dx 


ikM^X/0^ 


)]• 


(8) 


It  is  not  strictly  necessary  to  use  equations  (3)  and  (4)  to  calculate  the  under¬ 
lying  steady  flow  field.  An  alternative,  see  [4]  and  [5]  and  elsewhere,  is  to  use 
the  steady  solution  from  a  much  simpler  formulation  than  the  transonic  small 
perturbation  equation,  thus  meiking  computer  time  savings. 

Conversely,  the  steady  flow  field  could  be  generated  by  a  more  complicated 
equation  such  sis  the  full  potential  equation.  The  advantage  of  this  approach 
over  the  use  of  the  small  disturbance  equation  is  that  shock  strength  and 
position  are  modelled  more  realistically. 

An  attractive  feature  of  using  the  full-potential  equation  for  the  steady  flow 
field  is  that  this  equation  can  be  solved  with  relatively  little  extra  cost  over  the 
small  disturbance  equation.  Solving  the  steady  full  potential  equation  simply 
means  using  a  different  non-linesu-  term  <7,,  as  in  [15]. 
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The  method  of  solution  of  both  the  steady  and  unsteady  problems  is  btised  on 
the  Morino  integral  equation  method  [8],  in  which  Green’s  Formula  is  used  to 
construct  an  integral  equation  for  the  velocity  potential  (j>,  where  tj)  represents 
in  the  steady  case,  and  $u  in  the  unsteady  case. 

Green’s  formula  for  both  the  steady  and  unsteady  equations  relates  the  poten¬ 
tial  at  a  point  P  in  the  flow  domain  to  an  integral  of  the  potential  over  the 
aerofoil  boundary  C5,  the  wake  contour  Cw,  any  shock  discontinuities  Cj), 
and  the  contour  at  infinity  Coo-  It  has  the  form 

where  C  =  CsUCw^CdUCoo  (as  shown  in  Figure  1),  N  is  the  inward  normcJ 
to  the  fluid  in  the  Prandtl-Glauert  space,  and  a  represents  either  as  or  Ou. 

The  domain  function  E  is  defined  by 

{1,  ^p  within  flow  field 

5,  4>p  on  the  boundary  (10) 

0,  <Ap  outside  flow  field. 

For  the  steady  problem 

=  (11) 

(where  R  is  the  distance  to  the  singularity)  is  an  appropriate  Green’s  function, 
while  for  the  unsteady  problem 

G  =  ^iH^Q^hR)  (12) 

(where  Hq  '  is  the  Hankel  function  of  the  second  kind  of  order  zero)  is  an 
appropriate  Green’s  function. 

Using  the  far-field  expansions  in  [1]  and  [6],  tmd  the  Green’s  functions  above, 
it  can  be  shown  that  at  infinity  the  contribution  to  the  righthand  side  of  (9) 
is  at  most  a  constant.  Hence,  the  boundary  integral  over  Coo  may  be  omitted. 
Also,  the  integral  over  the  shock  discontinuity  Cq  can  be  combined  with  the 
volume  integral  in  such  a  way  that  there  is  no  explicit  contribution  to  <j>p  (as  a 
result  of  the  shock  jump  conditions)  [6],  [14].  Since  the  new  integred  equation 
generated  this  way  can  be  shown  to  differ  very  little  from  the  origin^ll  integr2J 
equation  “at  the  numerical  algorithm  level”  [14]  when  discretised,  the  original 
equation  can  be  used  ignoring  altogether  the  presence  of  shocks,  which  are  said 
to  be  “captured”.  As  in  [12],  artificial  viscosity  can  be  used  when  evaluating 
the  field  integrals  to  guarantee  the  absence  of  (i.e.  smooth  out)  shock  waves. 

A  primary  advantage  of  the  integral  equation  approach  is  that  no  complex, 
body-fitting  grid  generation  scheme  is  needed  within  the  flow  field  (such  as 
would  be  required  in  a  finite  difference  method).  Most  of  the  effort  in  grid 
generation  is  needed  at  the  aierofoil  boundary,  which  is  approximated  by  small 
elements  (line  segments  in  2D). 
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Although  the  flow  field  F  must  still  be  discretised  in  order  to  evaluate  the 
righthand  side  of  equation  (9),  a  cartesian  grid  (consisting  of  rectangular  pan¬ 
els)  is  sufficient  for  this  purpose,  and  this  grid  need  not  be  fitted  to  the  shape 
of  the  aerofoil  contours.  In  fact,  this  grid  need  only  be  generated  where  the 
transonic  terms  as  and  a^  are  significant,  i.e.  near  the  aerofoil  (see  [5], [15]). 
Also,  both  the  steady  and  unsteady  calculations  can  take  place  on  the  same 
grid. 


3.  PROPOSED  NUMERICAL  SCHEME  -  AN  OUTLINE 

The  numerical  scheme  proposed  here  involves  the  two  stages  described  above, 
namely,  initial  calculation  of  the  steady  flow  field,  followed  by  calculation  of 
the  unsteady  flow  field.  Each  stage  uses  the  same  grid,  but  with  differing 
Green’s  functions  and  a  terms.  The  calculation  scheme  for  each  stage  is  the 
same  and  proceeds  as  follows: 

Step  0  Evaluate  4>  =  4>C  on  the  aerofoil  contour  using  a  normal  (subsonic) 
panel  method,  that  is,  setting  the  field  source  strength  a  to  zero.  The  aerofoil 
surface  is  divided  into  panels  and  the  unknown  4>C  approximated  in  some 
fashion,  e.g.  as  being  constant  on  each  panel.  The  known  quantity  ^  is  also 
approximated  in  some  consistent  fashion.  The  equation  for  41Q, 

‘  (13) 

where  C  =  C5  U  Ctp  becomes  a  matrix  equation  for  the  discrete  4>q  values. 

Step  1  Evduate  <l>  =  <f>p  in  the  flow  field  F  using  the  values  of  calculated 
at  the  previous  step,  and  the  current  value  of  a.  The  flotv  field  is  discretised 
into  panels  on  which,  for  example,  <t>f  is  constant.  The  discrete  version  of 

yields  (bfSX  each  field  panel. 

Step  2  Evaluate  the  field  source  strength  a  by  applying  finite  differences  to 
formulas  (4)  or  (8).  The  expression  for  a  is  non-linear  for  the  steady  problem, 
while  for  the  unsteady  problem  it  is  linear  but  depends  on  the  calculated  steady 
solution. 

Exit  iteration  loop  if  field  source  strengths  have  converged. 

Step  3  Evaluate  4>  =  <t>c  on.  the  aerofoil  boundary  using  the  transonic  panel 
method  (with  field  source  term).  This  step  echoes  Step  0  and  involves  solution 
of  the  linear  matrix  equation  arising  from  discretisation  of 


Go  to  Step  1 
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Smoility  and  convergence  have  not  been  discussed  here  and  are  not  guaranteed. 
However,  in  [15]  encouraging  results  are  given  for  the  convergence  of  a  similar 
natural  iteration  scheme.  For  subcritical  (shock-free)  convergence  to 

acceptable  accuracy  occurs  within  5  iterations,  while  for  supercritical  flows 
(with  shocks)  convergence  occurs  within  30  iterations. 


4.  SUBSONIC  PROBLEM 


The  subsonic  problem  is  a  subset  of  the  transonic  problem,  corresponding  to 
Step  0  above,  in  which  the  field  source  strength  is  zero  and  the  steady  Green’s 
function  (11)  is  used.  This  problem  is  considered  first  in  order  to  identify 
areas  of  possible  difficulty  in  the  panelling  and/or  evaluation  of  the  influence 
coefficients,  and  to  provide  the  programming  framework  for  later  extensions  to 
steady  trtinsonic  flow  (involving  non-linear  terms  and  iteration)  and  unsteady 
flow  (more  complicated  Green’s  function). 

The  appropriate  boundary  value  problem  for  the  steady  potential  <{>  is  to  satisfy 
equation  (2)  with  associated  boundary  conditions 


9</> 

dn 


— i  •  n. 


(16) 


on  the  aerofoil  boundary.  As  seen  alread''.  the  Prandtl-Glauert  transformation 
puts  the  tranconic  small  perturbation  eqt.ation  into  the  form  given  in  equations 
(3)  and  (4),  while  the  boundary  conditions  (16)  become 


dp  _  i-N 

dN  ~  /?2  ■ 


(17) 


Note  that  N  is  the  normal  to  the  aerofoil  in  the  transformed  flow  domain, 
and  that  higher  order  terms  have  been  neglected  here,  consistent  with  the 
linearisation  of  the  governing  equation. 


The  integral  equation  (13)  for  p,  obtained  by  application  of  Green’s  formula, 
is  solved  by  the  method  of  collocation.  The  aerofoil  contour  Cg  is  assumed  to 
consist  of  Nt  discrete  panels  Pj  whose  end  points  are  (zj_i,  )  and  (xj,yj), 
with  equation 

x  =  ajy  +  hj.  (18) 

(In  [2]  the  wing  surface  is  fitted  exactly  and  the  boundary  integrals  are  carried 
out  numerically,  whereas  in  the  Morino  method  an  approximate  geometry  is 
used  and  the  integrtds  axe  usually  calculated  analytically,  especially  near  singu¬ 
larities.)  The  collocation  points  (fy,?;,)  are  taken  to  be  the  Nt  panel  midpoints, 
with  associated  values  of  the  potential  Pi  =  P(ii,rit). 


At  the  (th  collocation  point  we  have 

where  R,-  =  \/{x  -  +  (y  - 
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Assuming  o  is  constant  on  each  panel,  and  treating  the  wake  boundary  as  a 
straight  line  emanating  from  the  trailing  edge  across  which  there  is  no  pressure 
jump,  we  have 


dC 


(\ogRi)d( 


(20) 


for  i  —  1.2,...,  iV(,  where  (pte  Ihe  value  of  at  the  trailing  edge.  This  leads 
to 


=  E/> 


(21) 


Jztc 


R? 


for  i  =  1.2 . Kf.  where  (ite^yte)  the  coordinates  of  the  trailing  edge. 

and  Ntf  is  the  panel  number  of  the  (upper)  trailing  edge  panel.  Evaluating 
the  integrals  we  have,  finally 


St  St  A'(t  +  1 

=  Wijcpj 

j=l  j=l  }=S„ 

for  ;  =  1,2 . .V(,  where 

B,j  =  .2  -  t/i  +  aj(x  -  ^ilKlog  i?i  -  1) 


(22 


+  {ajrji  +  bj  —  ^,)arctan( 


y  -r]i  +  aj(x  - 


-4,7  = 


r  ^  ,y  -Vi  +  aj(x  -^i)  iV} 

-arctan( - ; - - - )'  for  a, 

[arct£in(-^^ — ^)]  ’ 

t  yi  ”  t/,-  Jrj_i 


for  a,  =  oo, 


and 


Wi 


'ij  =  ±  lim  arctan('^^^ — ^)  -  arctan(  — —)  ,  j  = 
A'— ool  yte-Vi  yie~mi 


x^te  ~  ^1 1 
Vie  -  rii ' 


(23) 


(24) 


(25) 


The  above  set  of  simultaneous  equations  is  now  solved  for  the  unknowns  tfij. 

In  the  unsteady  ca.'”,  the  Green’s  function  is  so  complicated  that  the  influence 
CO'  ificients  can  only  be  evtduated  numerically.  Numerical  experiments  have 
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determined  that  for  the  steady  problem  a  four  point  Gaussian  rule  is  sufficiently 
accurate  for  use  in  evaluation  of  the  influence  coefficients,  except  near  the 
collocation  point  where  it  becomes  necessary  to  extract  the  singularity  and 
evaluate  that  part  separately  (rmalytically ). 

For  interfering  aerofoils/components  the  wake  terms  W^j  need  very  careful 
treatment,  especicdly  when  the  idealised  (i.e.  straight)  wake  of  one  aerofoil 
component  intersects  with  another  component  downstream.  An  interpretation 
used  successfully  here  is  that  the  arctan  terms  in  M'lj  represent  the  angles  3 
and  a  respectively,  shown  in  Figure  2.  These  angles  change  continuously  as 
the  point  crosses  the  wake  (even  though  the  arctan  function  jumps). 

As  a  result.  ±27r  must  be  added  to  H',j  where  necessary  to  retain  continuity. 

Once  the  potential  has  been  obtained,  an  obvious  method  for  calculating  the 
pressure  coefficient  would  be  to  differentiate  equation  (22)  with  respect  to 
and  rji  and  use  these  derivatives  in  all  calculation  .  This  works  quite  well  within 
the  flow’  field  (provided  the  factor  of  tr  is  replaced  by  27r  in  (22)).  However,  it 
also  leads  to  large  errors  near  the  aerofoil  because  of  the  assumption  that  the 
potential  varies  as  a  step  function.  The  function  given  by  the  right  hand  side 
of  equation  (22)  becomes  tJmos,  step-like  between  collocation  points,  causing 
the  calculated  derivatives  to  be  useless  near  these  points,  i.e  on  the  aerofoil 
boundary. 

Since  direct  ctdculation  is  impossible,  an  alternative  is  to  calculate  derivatives 
by  finite  differences.  This  is  done  using  the  collocation  points,  at  which  0  is 
known  best. 


Returning  to  the  original  unsealed  co-ordinates,  the  pressure  coefficient  Cp  is 
calculated  from  the  following  formula,  applicable  to  thin  geometries; 

Cp  = -2[i  •  q-l- •  q]  4- .U^(i  •  q)^  (20 

w'here  q  =  Vo  is  the  perturbation  velocity.  On  the  aerofoil  boundary  we  write 


q  =  qcS  -I-  (jnii 


(27) 


where  s  and  n  are  the  tangential  and  normal  vectors  respectively.  The  pressure 
coefficient  can  thus  be  rewritten  as 


Cp  —  -  2(95(1  •  s)  +  -(95 


9n) 


-  2gsgi^(i  ■s)  +  q,  )], 


(28) 


which  has  the  advantage  that  the  magnitude  of  the  normal  velocity  component 
9n  =  -i  n  is  known  exactly  from  the  boundary  conditions,  while  the  magnitude 
of  the  tangential  velocity  9*  is  easily  approximated  by  finite  differences  at 
the  boundary.  Thus  finally,  the  following  formula  applies  for  the  pressiure 
coefficient  at  the  tth  collocation  point; 


,  ^  _  2  rA0;Ag,-  -I-  -  ^(At?,)^t 

^  i  (A^,•)2-|-(A9,■)•■^  J 

,  ,^2  -  2A,^,A^.(  A>7,)2  +  (At?;)'* 

^  [(A?.)2  +  (At,.)2]2 
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where 


A(?i  ~  (Pi+i  -  Oi-i, 

(30) 

=  'It-t-l  -  'li-l 

(31) 

Mi  =  ~  Cl-l- 

(32) 

5.  RESULTS  AND  DISCUSSION 

Results  have  been  computed  for  the  case  A/  =  0  for  three  test  configurations, 
namely 

1 .  A  circle  with  unit  radius  centered  at  the  origin. 

2.  A  XACA0012  aerofoil  at  zero  angle  of  attack. 

3.  Williams  configuration  B  [16]  (aerofoil  plus  flap  at  10  deg). 

Figures  3  to  5  are  plots  of  —  Cp  vs  chordwise  position  for  these  three  configura¬ 
tions.  The  exact  incompressible  solutions  are  denoted  by  lines  and  numerical 
solutions  are  denoted  by  asterisks. 

The  agreement  in  all  case-  is  excellent.  The  agreement  at  the  trailing  edges  in 
configuration  3  (Figure  5)  is  surprisingly  good  in  view  of  the  assumption  that 
each  wake  is  a  straight  line  emanating  from  the  trailing  edge. 


6.  CONCLUSION 

In  this  article  the  outline  of  a  numerical  scheme  for  the  calculation  of  two- 
dimensional  unsteady  transonic  aerodynamics  has  been  presented.  Excellent 
results  have  been  obtmned  for  the  first  step  towards  implementing  this  scheme, 
that  is,  for  the  special  case  of  steady  subsonic  aerodynamics. 

Despite  the  apparent  simplicity  of  the  outlined  scheme,  implementation  for 
unsteady  and  transonic  aerodynamics  will  still  involve  significant  leaps  in  diffi¬ 
culty.  For  unsteady  flow,  a  more  complicated  Green’s  function  will  mean  that 
the  integrals  involved  can  no  longer  be  evaluated  in  closed  form.  Difficulties 
are  also  anticipated  in  the  application  of  the  Kutta  condition  at  the  trailing 
edge.  For  tranronic  flow  the  addition  of  field  source  terms,  and  the  consequent 
need  for  8Ln  iterative  numerical  scheme,  will  lead  to  demands  on  computer  time 
and  storage.  Problems  with  stability  and  convergence  of  the  numerical  scheme 
are  also  a  possibility. 
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Figure  3. 
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